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Abstract 

We present an efficient and robust algorithm for finding points of 
collision between time-dependent parametric and implicit surfaces. 
The algorithm detects simultaneous collisions at multiple points of 
contact. When the regions of contact form curves or surfaces, it 
returns a finite set of points uniformly distributed over each contact 
region. 

Collisions can be computed for a very general class of surfaces: 
those for which inclusion functions can be constructed. Included 
in this set are the familiar kinds of surfaces and time behaviors 
encountered in computer graphics. 

We use a new interval approach for constrained minimization to 
detect collisions, and a tangency condition to reduce the dimension¬ 
ality of the search space. These approaches make interval methods 
practical for multi-point collisions between complex surfaces. An 
interval Newton method based on the solution of the interval lin¬ 
ear equation is used to speed convergence to the collision time and 
location. This method is more efficient than the Krawczyk-Moore 
iteration used previously in computer graphics. 

CR Categories: 1.3.5 [Computer Graphics]: Computational Ge¬ 
ometry and Object Modeling; G.4 [Mathematical Software]: Relia¬ 
bility and Robustness 

General Terms: collision detection, parametric surface, con¬ 
strained minimization, interval analysis 

Additional Key Words: inclusion function, interval Newton 
method, interval linear equation 

1 Introduction 

Detecting geometric collisions between curved, time-dependent 
(moving and deforming) objects is an important and difficult prob¬ 
lem in computer graphics. This paper discusses a practical and ro¬ 
bust algorithm for detecting collisions between objects represented 
as parametric or implicit surfaces. We ignore the problem of com¬ 
puting the physical response to collisions; much of this topic is 
treated in other work [BARA90,META92]. Instead, we concen¬ 
trate on the purely geometric problem of computing a solution set 
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Figure 1: Problem Statement: Given a collection of time-dependent curved 
surfaces, find a set of collision points representing the contact regions. In 
this example, the dots show the points detected by the collision algorithm 
when a torus moves down over a cone, contacting it in a circle. 

of points where a set of time-dependent surfaces first contact (Fig¬ 
ure 1). 

Previous work on geometric collision detection is fairly extensive, 
both in computer graphics and in other fields such as CAD/CAM 
and robotics. Detection of collisions between polyhedral objects 
was studied in [MOOR88]. Baraff [BARA90] presented a method 
of computing collisions between parametric or implicit surfaces 
by computing extremal points using non-linear equation solvers. 
Sclaroff and Pentland [SCLA91] present a method for detecting 
collisions between implicit surfaces by “plugging” vertices of a 
polyhedral approximation of one surface into the inside-outside 
function of the other. Von Herzen, et. al., [VONH90] presented an 
algorithm for detecting collisions of parametric surfaces using Lip- 
schitz bounds. Duff [DUFF92] used interval methods to compute 
collisions between boolean combinations of implicit surfaces. 

To make collision detection practical, much of the previous work 
traded off accuracy and robustness for efficiency, or limited the 
kinds of shapes that could be handled. Polyhedral methods such 
as in [MOOR88], although fairly efficient, are not well suited to 
surfaces that deform in time. Exploiting coherence for rolling or 
sliding contact of polyhedral objects is difficult, and use of a fixed 
sampling mesh can cause severe approximation errors. Polyhedral 
methods also require many numerically difficult special cases which 
led [MOOR88] and [SCLA91] to neglect cases where “tunneling” 
may occur either between polygon edges or between small implicit 
surfaces passing entirely through a large polygon. 

Baraff [BARA90] chose to limit objects to the union of con- 
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take one or more steps in the ODE solver 
compute collisions in the resulting time interval 
if a collision occurs (at time f) in the interval 
compute a collision response 
reset ODE solver to t = t* 

endif 


Figure 2: Computational Model for Collision Detection and Response 

vex polyhedra and strictly convex closed surfaces. This restriction 
simplified his collision detection algorithm and allowed tracking 
of single contact points between curved objects. He did not treat 
non-convex surfaces (such as saddle shapes) and manifolds with 
boundary (such as half a sphere). We solve the problem for a more 
general class of surfaces with many points of contact, as shown in 
Figure 1. 

As noted in [VONH90], methods which depend solely on point- 
wise evaluations, including the above methods, cannot guarantee 
accurate collision detection. To solve this problem, Von Herzen 
bounded the output of functions over a region using a Lipschitz 
bound. Duff [DUFF92] used interval analysis to produce tighter 
bounds than Von Herzen’s Lipschitz bound. Both of these methods 
used binary subdivision to search for collisions; we speed up the 
approach significantly by combining binary subdivision with an 
interval Newton method. 

The technique we describe offers several fundamental improve¬ 
ments over previous techniques: 

1. The most novel aspect of our technique is the ability to detect 
simultaneous collisions (multiple contacts at the same time), 
even when the collisions occur at a higher dimensional mani¬ 
fold of contact, rather than at a set of isolated points. In this 
case, the algorithm samples the region of contact with a finite, 
uniformly-distributed set of points. The spatial sampling den¬ 
sity is a parameter to the algorithm. To our knowledge, no 
previous algorithm handles this situation. 

2. Our technique works for both rigid and deforming objects, and 
for implicit or parametric objects. 

3. Our technique is practical for computer graphics applications, 
and has been used in animations involving hundreds of objects. 

4. Our technique includes a method (tangency constraints) to 
reduce the dimensionality of the space of possible solution 
points, as shown in Figure 3, dramatically speeding up the 
method. The tangency constraints also provide a square system 
of equations for the interval Newton method, helping us detect 
isolated point collisions. 

5. Our technique uses a test for uniqueness of roots of a system of 
equations in a region. This test can be verified in many cases, 
allowing the algorithm to terminate without further subdivision 
around collision points. 

6. Our technique can be used both to compute collisions between 
formerly disjoint bodies which come into contact, or to com¬ 
pute additional points of contact between bodies as they roll or 
slide over each other (see Section 1.1). 

1.1 Fitting Collision Detection into a Larger System 

Figure 2 shows how collision detection fits into a larger pro¬ 
gram for computing physical simulations of dynamic systems. The 
system is composed of three parts: the ODE (ordinary differential 
equation) solver module, the collision detection module, and the 


collision response computation module. The ODE solver computes 
the motions of objects over time, using equations governing the dy¬ 
namic behavior of bodies, and produces a functional representation 
of the motion. 1 Motion is computed without considering collisions, 
so that the results are only valid until the next collision occurs. The 
collision detection module takes the functional representation pro¬ 
duced by the ODE solver and computes when and where the first 
collision occurs in the given time interval. If a collision occurs, a 
collision response is computed, which may discontinuously change 
the state of the system of bodies. The ODE solver continues forward 
in time from this computed collision time, discarding any state after 
it. 

Two modes of operation are required in collision detection: 

1. compute any collisions for bodies that are initially not in con¬ 
tact 

2. compute additional collisions for bodies that are already in 
continuous (rolling or sliding) contact 

The algorithm described in this paper handles both situations. For 
greatest efficiency and modularity, we advocate handling coherence 
in the ODE solver. By coherence, we mean the tracking of contact 
points between bodies rolling or sliding over each other. In these 
situations, collision detection is required only to compute new points 
of contact not already tracked by the ODE solver (mode 2 above). 
The solver must therefore inform the collision detection module of 
the motion of the contact points it is tracking, so that these points 
may be excluded from consideration (see Eq. 7). The collision 
detection module must also compute the initial points of contact 
when the simulation is begun or when continuous contact begins 
between bodies (mode 1 above). 

1.2 Overview 

The mathematics of the collision detection problem is treated in 
Section 2. Sections 3, 4, and 5 discuss the constrained minimiza¬ 
tion algorithm, an interval Newton enhancement, and termination 
criteria, respectively. Section 6 presents a simple culling test which 
discards non-colliding surface pairs and tightens a bound on the 
collision time. The full collision algorithm, combining constrained 
minimization, the culling test, and other tools from computational 
geometry, is presented in Section 7. Our technique, like all interval 
methods, requires inclusion functions, whose construction is sum¬ 
marized in Section 8. Finally, results and conclusions are described 
in Sections 9 and 10. Appendix A extends our approach to sur¬ 
faces that are piecewise smooth by adding conditions for face, edge, 
and vertex interactions (see Figure 11). Appendix B describes the 
construction of inclusion functions for Chebyshev polynomials. 

2 The Collision Problem 

The equations that specify that two surfaces collide may be di¬ 
vided into two parts: a contact constraint, that specifies that the 
two surfaces intersect, and a tangency constraint, that specifies that 
the two surfaces are tangent at their point of intersection. The tan¬ 
gency constraint reduces the dimensionality of the space of possible 
collision points, as shown in Figure 3. It also allows faster conver¬ 
gence (using interval Newton, which we will describe in Section 4) 


1 In our rigid body simulations, the solver produces a time-varying quaternion and 
translation vector. Each component of the quaternion and vector is represented using 
univariate Chebyshev polynomials. 
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is the zero vector): 



contact with tangency contact without tangency 


Figure 3: Reducing the Dimensionality of the Space of Collision Points 
Using the Tangency Condition: The intersection of two bodies (like a sphere 
moving to the right with a stationary plane) typically forms a whole 2D 
manifold of contact through time. With the tangency constraint, the solution 
space is often reduced to one or a few points by eliminating cases like that 
shown on the right. Reducing the solution space to an isolated space-time 
point is one of the ideas that makes this method practical. 



incoming collision 



outgoing collision 


Figure 4: Incoming and Outgoing Collisions: The unbroken circles repre¬ 
sent bodies later in time. A dot represents the collision point; the arrows 
represent the direction of movement. 

and robust testing of isolated collisions (using an interval solution 
uniqueness test described in Section 5). 

We also distinguish between incoming collisions, in which the 
surfaces collide by moving closer to each other, and outgoing colli¬ 
sions, in which the surfaces are interpenetrating and become tangent 
as they move apart. These situations are compared in Figure 4. The 
distinction is necessary in the simulation of dynamic systems where 
each surface encloses a solid. Eliminating outgoing collisions al¬ 
lows the simulator to ignore collisions which were previously de¬ 
tected; i.e., collisions between surfaces already in contact which are 
moving away as a response to the collision. 

2.1 Parametric Surfaces 

Let two deforming parametric surfaces be represented by the 
twice-differentiable mappings S\(u\, vi, t) and S 2 (u 2 , V 2 , t), where 
Sf. R 3 —► R 3 . At a particular instant of time, each of the surfaces 
is formed by the image of Si over a rectangle in (ui,vi) space. 2 
In this section, we consider the case of collisions between solids 
each bounded by a single, smooth, closed parametric surface. Ap¬ 
pendix A generalizes the discussion to parametric surfaces which 
are only piecewise smooth. 

Contact Constraint The contact constraint merely states that the 
two surfaces intersect (i.e., the vector difference of the two surfaces 


2 Using a rectangular domain for parametric surfaces does not limit the kinds of 
surfaces that can be collided. Parametric surfaces defined on non-rectangular domains 
can be handled by mapping a rectangle into the required non-rectangular domain before 
mapping onto the surface [SNYD92b]. 


Sl(u\, Vl, t) — £2(^2, V2, 0 = 0. (1) 


Tangency Constraint The tangency constraint implies that the 
instantaneous normal vectors on the two surfaces at their point of 
contact are anti-parallel. Stated another way, the (u, v) tangent 
vectors on one surface must be perpendicular to the instantaneous 
normal vector on the other surface. We thus have the following 
system of two equations 3 


^ Vl, t) • N 2 (U 2 ,V2, t) 

OU\ 

. dSi 

\ —— - N 2 (u 2 ,V2, t) 

\ OV 1 * 


( 2 ) 


where N\ and N 2 are the outward normal vectors to the surfaces S\ 
and S 2 , respectively, given by 


Niim, Vi, t) = —^{ui, Vi, t) x —^(ui, Vi, t) for i- 1,2. 

uUi (JVi 

The algorithms that follow here assume that N\ and N 2 are nowhere 
0; that is, surfaces have a nonvanishing normal vector everywhere 
and for all relevant time. 4 The whole collision equality constraint 
is given by a nonlinear system of 5 equations in 5 variables, three 
from Eq. 1 and two from Eq. 2. 


Incoming Constraint The incoming collision condition states 
that the relative velocity of the collision point must face the same way 
as the surface normal (the two vectors must form an acute angle), 5 
and the two normals must face in opposite directions (forming an 
obtuse angle). This condition yields two inequality constraints: 


(^P~(ui,vi,t) - v 2 , 0 ) -N\{ui,vi,t) > 0 

dt dt (3 ) 

and — N\(u\, vi, t) • N 2 (u 2 , v 2 , t) > 0. 


2.1.1 Example: Rigid Parametric Surfaces 

The above constraints may be applied to the special case of rigid 
parametric surfaces. In this case, we have two time-independent 
surfaces s\(u\, vi) and ^ 2 (^ 1 , vi). The time-varying version of these 
surfaces is given by 

Si(ui, Vi, t ) = Ri(t) Si(u\,v\) + Ti(t ) for i = 1,2 

where Ri(t) is a time-varying rotation matrix and Ti(t) is a time- 
varying translation vector, specifying the trajectory of surface V s 
coordinate origin. 

Contact Constraint The contact constraint may be expressed as 
R\{t) s\{u\,v\) + T\(t) R 2 (t) s 2 (u 2 , v 2 ) - T 2 (t) =0. 


3 A similar, though functionally dependent, constraint may be derived by switching 
Si and S 2 . 

4 If the calculated normal vector becomes zero, such as at the poles of a parametric 
sphere, the tangency constraint becomes trivially true. The algorithm will therefore 
rely on the contact constraint to detect a collision in this case. 

5 We assume here that the surfaces are parameterized so that the normals N\ and N 2 
face outward. 
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Tangency Constraint Let m(u\, vi) and m(u 2 , v 2 ) be the time- 
independent normals of the surfaces s\ and 52 , given by 

, X V V 

rii(ui,Vi ) = —(m/, V/) x —(ui,Vi). 

C/Ui t/V/ 

The time-varying surface normals can therefore be expressed as 


#i(Mi, v/, 0 = i?i(0 rii(ui , Vi) 


since i?/0 is a rotation matrix. The tangency constraint is then given 

by 

_ q 

^( R l(t) ■ N 2 (u 2 ,V2,t) J 

Incoming Constraint The incoming constraint is given by 

[ki(t)si(u u Vl) + Tl(0 - R 2 (t)s 2 (u 2 , v 2 )— 

T 2 (0] • M(mi, vi, 0 > 0 

and — N\(ui,vi,t) • N 2 (u 2 ,V 2 ,t) > 0 

where and 7) are the time derivatives of the rotation matrix and 
translation vector of the two surfaces. 


Figure 5: Simultaneous Collisions — The tori collide at two isolated points. 

We also note that CSG operations on implicit surfaces, as in 
[DUFF92], can be handled very efficiently using our techniques. 
Assume the surface F\ is represented as the boolean subtraction of 
two simple implicit surfaces F a (x, y, z, t) and F)?(x, y, z, t). The first 
equation in the contact constraint (Eq. 4) then becomes 

(F a = 0 and 7/ > 0) or (Ft = 0 and F a < 0) 



2.2 Implicit Surfaces 

Let two time-varying implicit surfaces be represented using the 
scalar functions F i(x, y, z, t) andF 2 (x, y, z, t). Points on each surface 
are defined as the zero-sets of these functions. 


Contact Constraint 

equations 


The contact constraint is the system of two 


(Fi(x,y,z,t)\ „ 
\F 2 (x,y,z,t) ) 


(4) 


Tangency Constraint Let the function V7q(x, y, z, 0 be the spatial 
gradient of the implicit functions (i.e., with respect tox, y, and z). 
The tangency constraint is then given by 


assuming implicit surface functions are positive outside the surface 
they represent. Similar restricted equality constraints can be derived 
for the tangency constraint. 

Constraints for rigid motion of implicit surfaces are easily de¬ 
rived by applying the above general equations to the rigidly moving 
implicit surface 

F(x, y, z, t) = f(w) where w = R r (t)((x, y, zf - T(t )) 

where/: R 3 —► R is the implicit equation of the time-independent 
surface, R T (t) is the transpose of the time-varying rotation matrix, 
and T(t) is the time-varying translation vector. 

2.3 Collision As a Constrained Minimization Problem 


V F\ (x. y, z, t) x VF 2 (x, y, z, t) = 0. (5) 

This constraint, although a system of three equations, contains only 
two functionally dependent equations. The entire collision equality 
constraint for implicit surfaces is thus given by a system of five (four 
functionally independent) equations in the four variables x, y, z, and 
t (Eqs. 4 and 5). 6 

Incoming Constraint The incoming constraint is given by 


--^-(x, y ,z,t)\\VF 2 (x,y,z,t)\\ - 

~^~(x,y,z, t) \\VF x (x,y,z, Oil > 0 ( 6 ) 

and - VFi(x, y, z, t) • VF 2 (x, y, z, t) > 0. 


6 We note that detecting collisions between implicit and parametric surfaces is a 
simpler problem than colliding pairs of parametric or implicit surfaces. By substituting 
the output of the parametric surface as the input ( x, y, z) of the implicit surface, a system 
in 3 variables, (u, v, t ), results, where u and v are the parametric surface coordinates. 


The final collision constraint may be described as an equality con¬ 
straint involving a function C (for the contact and tangency con¬ 
straints) and a logical composition of inequality constraints in¬ 
volving a function D (for the incoming constraint). For colli¬ 
sions between parametric surfaces, C and D are vector functions 
of (u\, vi, w 2 , v 2 , t) (equations 1-3); for implicit surfaces they are 
functions of (x, y, z, t) (equations 4-6). We are interested only in the 
minimum t collision, since our representation for the time behavior 
of the surfaces may be invalid after this time. 

The desired collision time for parametric surfaces, t *, can there¬ 
fore be expressed using the constrained minimization problem 


minimum 

Ol ,vi,U 2 ,v 2 ,t)eX 0 



C(u\ , Vl, U 2 , v 2 , 0 = 0 
D(u\ , vi , w 2 , v 2 , 0 > 0 


and 


A similar statement results for detection of collisions between im¬ 
plicit surfaces. We would like to compute t* or detect that the 
constraint is satisfied nowhere in the parameter space Xq. 

We also need the location of the collision and the surface normal 
vectors there. There may be multiple points of contact at the time 
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of collision, which we call simultaneous collisions, as shown in 
Figure 5. The points of contact may be a finite number of isolated 
points as in Figure 5, or they may form a curve or surface, called the 
contact manifold. For example, if the falling torus from Figure 5 
were in the same orientation as the stationary one at the bottom, the 
collision points would form a circle. 

To detect simultaneous collisions, we need to detect minimum 
t solutions which are simultaneous or simultaneous within some 
tolerance c. Mathematically, we require the set of collision points 
(u \, v*, u 2 , vj) such that 

Cyui ,Vi,u 2 ,v 2 ,t ) = 0 and D{ui ,Vi,u 2 ,v 2 ,t ) > 0. 

For a given collision point, the location of the collision, p*, is 

p* = Si(u *, v* , t*) = S 2 (u 2 , v 2 , t*). 

The normal vectors at the collision may be defined similarly by 
evaluating N\ and N 2 at the collision point. 

To compute a collision among a set of N time-varying parametric 
surfaces Si(ut, v/, t), we first use simple culling procedures to exclude 
pairs of surfaces which can’t collide. Section 3.2 will discuss a 
method of solving these sets of constrained minimization problems 
which can compute simultaneous collisions and which does not 
spend undue computation on collisions which occur after f . It 
returns the collision points when these occur at a finite set of isolated 
points, or a finite subset of the collision points uniformly distributed 
over the contact manifold. 

3 Interval Tools for Computing Collisions 

We now turn to a discussion of the interval tools necessary to solve 
the sets of constrained minimization problems that arise in colli¬ 
sions. 

3.1 Review of Interval Analysis 

An interval , A = [a, b], is a closed subset of R defined as 

[a, b] = {x | a < x < b, x, a, b E R} . 

The lower and upper bounds of an interval are written as 

lb \a,b\ = a 

ub \a,b] = b. 

A vector-valued interval of dimension n,A = (Ai, A 2 , ... ,A n ), is a 
subset of R n defined as 

A = {x\ Xi e Ai,i= 1,2,...,«} 

where each A/ is an interval. An interval A/ that is a component of 
a vector-valued interval is called a coordinate interval of A . 

The width of an interval, written w([a, b]), is defined by 

w([a, b]) = b — a. 

The midpoint of an interval, written mid([a, b ]), is defined by 

mid([a, b]) = Af- 

Similarly, the width and midpoint of a vector-valued interval of 


dimension n , A, are defined as 

w(A) = max w(A ; ) 

/=1 

mid(A) = (mid(Ai), mid(A2),..., mid(A n )). 

Hereafter, we will use the term interval to refer to both intervals 
and vector-valued intervals; the distinction will be clear from the 
context. 

An inclusion function for a function /, written □/, produces an 
interval bound on the output of / over an interval representing its 
input domain. Mathematically, for all intervals X in the domain of 
/, if a point x is in the input interval X then fix) is contained in the 
output interval □/(X); i.e., 

x E X => fix) E D/(X) for allx E X. 

Much more information about inclusion functions and their prop¬ 
erties can be found in the literature (see, for example, [MOOR79, 
ALEF83,RATS88]). Section 8 and the Appendices discuss ways to 
create inclusion functions given the functions they are to bound. 

3.2 Constrained Minimization Algorithm 

The constrained minimization problem involves finding the global 
minimizers 7 of an objective function/: R" R for all points that 
satisfy a constraint function F: R n —► {0,1}. 

For the case of computing collisions between parametric surfaces, 
we have the following variables, objective function, and constraint 
function: 

X = (mi, Vi, U2, V2, t) 

fix) = t 

F{x) = (C(wi, vi, U2, V 2 , t) = 0) and 
(D(u\, VI, U2 , V2, t) > 0 ) 

A region in the minimization algorithm is a 5D interval vector of 
the form 


X = 

(U u Vi,U 2 ,V2,T) 

= 

( [«UUvUH«UUvUH^1) 

where the superscripts l and u denote lower and upper bounds. The 
relevant inclusion functions are 8 

q/oo = 

VA 

□ F(X) = 

r [0,1], if lbQC(X) < 0, ub QC(X) > 0, 
l and u b QD{X) > 0 

( [0,0], otherwise 


where dC is an inclusion function for the collision equality con¬ 
straint C, and DD is an inclusion function for the incoming inequal¬ 
ity constraint D. 

The algorithm in Figure 6 finds solutions to the constrained min¬ 
imization problem in a specified region Xo. The algorithm uses a 
priority queue to order regions based on the upper bound of the 
objective function. Regions bounding the set of global minimizers 


7 The global minimizers are the domain points at which the global minimum of the 
objective function is achieved, subject to the constraints. 

8 The terminology lb □C(X) < 0 denotes that lb □ Q (v) < 0 for each component 
interval i of UC{X) (and similarly for upper bounds). 
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Minimize(a/,aF,A,X 0 ,a</e/) 
place initial regionX 0 on priority queue L 
initialize fs upper bound u <— +oo 
initialize solution set S <— 0 
initialize singular solution set S <— 0 
while L is nonempty 

get next region F from L 
if lb □/(F) > u + ( discard Y 
else if lb □/(F) > u - e and 

there exists N e S such that || od(Y) - □c/(5,)|| < 8 
then discard F 

else if F satisfies acceptance criteria A then 
add F to solution lists 
if F doesn’t contain a unique feasible point 
add F to S 

endif 

u <— min(w, ub C/(F)) 

delete from S and S all S t 9 lb □/(#) > w + e 

else 

subdivide F into regions Fi and F 2 
for F/ g {F,.F 2 } 

evaluate of on F; 
if nF(F/) = [0,0] discard F* 
evaluate □/ on F, 
if lb □/(/•) > u + e discard F* 
insert F* into L according to ub □/(/•) 
endfor 

endif 

endwhile 

Figure 6: Global Constrained Minimization Algorithm: This algorithm 
finds the global minimizers of an objective function /, with constraints F, 
acceptance criteria A, initial region Xo, solution distance mapping function 
d, simultaneity threshold e, and solution separation distance 8. 

are subdivided until they are rejected or satisfy the acceptance cri¬ 
teria, A, and are accepted as solutions. It halts with an empty list of 
solutions if there are no solutions to the constraint function in Xo, 
or a list of regions, S , representing the set of global minimizers of 
the constrained minimization problem. 

The variable u is a progressively refined least upper bound for the 
global minimum of the objective function. If we were only looking 
for a single collision point, we could halt the algorithm immediately 
after finding the first solution. To find collisions at multiple points 
of contact, the algorithm must be continued until the priority queue 
is empty. The variable u helps to prune the search after finding the 
first solutions. 

Selecting Finite Sets of Points from Contact Manifolds The 

parameters c, 8, and □ d allow the algorithm to select a finite set of 
regions distributed “uniformly” within the set of global minimizers, 
when this set is not finite. The parameter e is the simultaneity 
threshold , which specifies how close the value of the objective 
function must be for two points to be considered global minimizers. 
For collisions, t specifies how close in time two events must be in 
order to be considered simultaneous. The parameter 8 is the solution 
separation distance , which specifies how far apart two accepted 
regions must be to be accepted as separate solutions. The parameter 
□ J is an inclusion function for the mapping which takes points in 
parameter space to points in whatever space we desire distances to 
be compared. We call the function d the solution distance mapping 
function. 

As the algorithm progresses, it maintains two solution lists, S and 
S. S contains all accepted regions. We call S the singular solution 
set. The elements of S not in S are regions in which the existence 
of a unique feasible point has been verified. The statements 


if lb □/(F) > u - c and there exists Si e S 
such that || nd(Y) - □rf(S / )|| < 8 
then discard F 

check that the region F is not too close to regions already accu¬ 
mulated onto S. Note that the test I b D/(F) > u — c is critical to 
ensure that F doesn’t have an objective function value small enough 
to invalidate all the currently accepted regions. 9 

We use two lists, S and S , so that in the case that the global 
minimizers form a finite set of points, the algorithm can find all 
such points without discarding some based on distance to those 
already found. The algorithm is therefore able to resolve multiple 
isolated collisions that happen in a small area, regardless of the 
value of 6. 10 

Ordering Based on Upper Bounds Constrained minimization 
algorithms that have appeared before [RATS88,SNYD92c] order 
regions based on the lower hound of the objective function. We use 
the upper bound to make tractable computing solutions on a contact 
manifold. 

At any time, the union of all regions on the priority queue forms a 
bound on the set of global minimizers of the constrained minimiza¬ 
tion problem. As the algorithm progresses, regions are subdivided 
or rejected, so that the regions which remain on the priority queue 
become a tighter bound on this set. Because of the inclusion mono¬ 
tonicity property of inclusion functions, 11 as regions on the queue 
shrink, the computed lower bound on the objective function tends 
to increase and the upper bound tends to decrease. 

Assume the set of global minimizers forms a continuous manifold 
rather than a finite collection of isolated points, as shown Figure 1. 
If the priority queue is ordered using lower bounds, when a given 
region is subdivided, its children will generally have larger lower 
bounds for the objective function, and will be placed in the priority 
queue behind less highly subdivided regions. A breadth first traver¬ 
sal tends to result, with less highly subdivided regions examined 
first. If we have a whole manifold of global minimizers and strin¬ 
gent acceptance criteria, we will have to compute a huge number of 
tiny regions bounding the entire solution manifold before even the 
first region is accepted as a solution. 

By ordering based on the upper bound, more highly subdivided 
regions tend to be examined first because they tend to have smaller 
upper bounds. We quickly get to a region which is small enough to 
satisfy the acceptance criteria. This allows our upper bound u to be 
updated. It also allows regions to be accumulated onto our singular 
list S. Regions that are too close to any member of S can then be 
eliminated, making it possible to find a distribution of points on the 
contact manifold without undue computation. 

Acceptance Criteria The constraint inclusion function, DF(X), 
because it contains an equality constraint, returns either [0, 0] (i.e., 
the constraint is satisfied nowhere in X) or [0,1] (i.e., the constraint 


9 If a region can possibly have a feasible point with a value of / less than e from the 
value off in regions on S, we should not reject it just because it is close with respect to 
the function d to these regions. The algorithm might then discard a global minimi z er 
because of its closeness to regions which are possibly far from the global minimizer, in 
terms of bounds on the objective function. 

10 One problem with this technique is that if the collisions happened at a contact 
manifold and a finite number of additional isolated points, the algorithm may discard 
some of the isolated points because of the closeness criterion. We consider this problem 
minor since the set of global minimizers is infinite and the algorithm must chose a subset 
anyway. 

11 An inclusion function, □/, is inclusion monotonic if F C X □/( Y) C □/(X). 
In practice, the standard ways of constructing inclusion functions generate inclusion 
monotonic inclusion functions. 
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may be satisfied in X). We must resort to other means to determine 
if the constraint is actually satisfied. Section 5 discusses conditions 
which guarantee that a region contains a unique solution to the 
equality constraints. These conditions can therefore be used as 
acceptance criteria in the algorithm, which we call the isolated 
point acceptance criteria. They also allow the upper bound u to be 
updated via 

u <— min(w, ub □/(F)) 

since Y is guaranteed to contain a feasible point. 

The algorithm also makes use of emergency acceptance criteria 
which do not guarantee a unique solution but are guaranteed to 
be satisfied for regions of small enough size. 12 The simplest such 
criterion is w(X) < e; a better one is w( □£/£)) < e where OS is the 
parametric mapping of one of the colliding surfaces. Regions which 
are accepted via the emergency acceptance criteria are inserted both 
onto the list of solutions, S , and the singular solutions, S. 

Subdivision The simplest method of subdividing candidate in¬ 
tervals in the minimization algorithm is bisection, in which two 
intervals are created by subdividing one of the input dimensions at 
its midpoint. Many methods can be used to select which dimension 
to subdivide. For example, we can simply pick the dimension of 
greatest width. A better alternative is to scale the parametric width 
by some measure of its importance to the problem we are solving. 
For each variable in the collision problem x,-, and a given candidate 
region X, we have used a scaling value st defined by 

Si= y^max(|lbn||(X)|,|ubn||(X)||). 
j= 1 

Here, / refers to the equality constraints (contact and tangency) of 
the collision problem. We then pick a dimension to subdivide, i, 
such that the scaled width Si w (X/) is largest. 

Given a candidate interval, techniques also exist which allow 
us to compute a smaller interval which can possibly contain feasi¬ 
ble points of the constraint. These methods and how they can be 
added to our simple minimization algorithm are discussed in Sec¬ 
tion 4. Even more sophisticated subdivision methods exist, such as 
Hansen’s method which involves accumulating gaps inside candi¬ 
date intervals by using infinite interval division (see [RATS88] for 
a full description). 

Multiple Element Constrained Minimization The algorithm 
can easily be modified to accept an array of sets of minimization pa¬ 
rameters (□/, OF,A,Xo)i. This allows simultaneous solution of sets 
of problems from different pairs of surfaces, or different tangency 
situations for the same pair of piecewise parametric surfaces. As a 
result, computation is not wasted on collisions which happen after 
the first collision, t > t* . We call this modified constrained mini¬ 
mization algorithm the multiple element constrained minimization 
algorithm. 

Sets of minimization subproblems may be implemented by asso¬ 
ciating the array index of the appropriate minimization subproblem 
with each region inserted onto the priority queue, and using the ap¬ 
propriate indexed inclusion functions and acceptance criteria when 
processing the region. 


12 Although we cannot guarantee a region X contains a solution, we can guarantee 
that it is arbitrarily close, in the sense that w(DC(Z)) < e where DC is an inclusion 
function for the collision equality constraint function C. 


Avoiding Detection of Tracked Points We can add additional 
inequality constraints to the constraint function F in order to avoid 
detecting collisions which occur at contact points already being 
tracked. If p is such a tracked point on a surface S(u, v, t ), the ODE 
solver computes a trajectory for p - S(u(t ), v(t), t). We then dis¬ 
card all global minimizers to the constrained minimization problem 
which satisfy 

\\S(u, v, t) - S(u(t), v(t), 0|| < A, (7) 

where A is a constant chosen by the user. The functions u and v 
have known representations, as computed by the solver. A natural 
interval extension of this constraint involving an inclusion function 
for S is then included in the constraint inclusion OF. An additional 
constraint is added for each tracked point. 

4 Interval Newton Methods 

In order to more quickly refine our intervals towards the solutions 
of the collision equality constraint C = 0, we make use of an 
interval Newton method. Interval Newton methods are applicable 
to the general problem of finding zeroes of a differentiable function 
/: R n R m in an interval X C R”. They allow us to find an interval 
bound on the set 

X* = {x E X | f{x) = 0} 

Let Z(X) be such a bound (i.e., X* C Z(X)). We can reduce the size 
of our candidate region X by 13 

x'=x pi z(x) 

In particular, Z(X) f] X = 0 implies that X contains no solutions. 
We call the operator Z(X) the interval Newton operator. 

Since X can only decrease in size after it is intersected with Z(X), 
this procedure can be applied iteratively to produce smaller and 
smaller regions, as in 

Xi +l =(Xip\Z(Xi)) 

Note however that a smaller region is not necessarily produced. 
Interval Newton methods should therefore be combined with bisec¬ 
tion. When interval Newton iteration is effective at reducing the 
size of X its use is continued. Otherwise, bisection subdivision is 
performed. 

The following sections present three methods for computing 
Z(X): 

• use of the Krawczyk-Moore form (Section 4.1) 

• use of the interval inverse (Section 4.2.1) 

• use of matrix iteration (Section 4.2.2) 

In each case, we modify the constrained minimization algorithm 
from Figure 6 by replacing the subdivide step with the code shown 
in Figure 7. 

4.1 Fixed Point Methods: the Krawczyk-Moore Form 

The familiar (point-wise) Newton’s method is used to converge on 
the solution to a system of equations/(x) = 0 where/: R n —► R n . 


13 A Pj B denotes the interval formed by the intersection of the intervals A and B. 
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Newton(y) 


compute interval Newton step ony 
if step succeeds then 

y' <- yp|Z(y) 


if Y' = 0, discard Y 

(proceed with next region) 

else if y' is sufficiently smaller than Y, insert Y 1 into L 
(proceed with next region) 
else subdivide Y ! 

(continue with Y replaced by Y f ) 
else subdivide Y 

(continue with Y) 


Figure 7: Interval Newton Modification to the Constrained Minimization 
Algorithm: The above algorithm replaces the subdivide step in the algo¬ 
rithm of Figure 6. 


The method starts with an initial “guess” at a solution xo and iterates 
via 

Xi+ i = pto) 


where 

p(x) =x- Yf(x). 


Y is a nonsingular n x n matrix which in straightforward Newton’s 
method is the inverse of the Jacobian matrix of / at x, i.e. 


Y = J~\x) 


Under certain conditions, this iterative procedure converges to a 
fixed point x*. If convergence is achieved, then the fixed point x* is 
a solution since 

pix*) — x* <=^>/(x*) = 0 
because Y is nonsingular. 

An interval analog of this method may be developed. Let X be an 
interval in R" in which zeroes of / are sought. We require a bound 
onX*. But 

V = {xex\f(x) = 0} 

= {xex \ P (x)=x} c (□/>(*)p|X) 

where □ p(X) is an inclusion function for the Newton operator p(x). 

The Krawczyk-Moore form, K(X, c, y), provides the necessary 
inclusion function for the Newton operator pix). It is simply a mean 
value form for p (see [SNYD92c] for a discussion of the mean value 
form) given by 


4.2 Linear Interval Equation Methods 

A second method for finding an interval bound on X* involves 
solving the interval analog of a linear equation. 

Let the coordinates of x be xi, X 2 ,..., x n . By the Mean Value 
Theorem, given a c G X, for each x G X, there exist n points, 
. £2 .such that 

fix) =/(c)+/(6,...,6)(x-c), 

where the jacobian matrix J is given by 

Jijdfi2 , . . . ,-;§k) = TfklO 

dXj 

and where each £* G X. Let □/ be an inclusion function for the 
Jacobian matrix of/, i.e., 

aj(x) = S^j \ Jij e D^.(X) 

If x is a zero of/, then there exists J G n/(X) such that 

fix) - 0 =/(c) + J(x — c). 

Therefore, if Q(X) is the set of solutions 

Q(X) = {x | f(c) +J(x-c) = 0 for some J G □/(X)} , 

then Q{X ) contains all zeroes of / in X. 

To compute an interval bound, Z, on Q(X), let y = x — c, and let 
Z' be an interval bound on the set 

{y \ Jy = —fie) for some J G D /(Y)} . 

Then the interval Z defined using interval addition as 

Z = Z ; + [c, c], 

is an interval bound on QiX). Thus, computing the interval Newton 
bound Z can be accomplished by solving an interval linear equation 
of the form 

Mx = b 

where M = □ J{X) is an n x n interval matrix, and b = —/(c) is an 
interval vector. 14 Stated another way, we require a bound on the set 

QiM, b) = {x | 3M G M. g b such that Mx . 


KiX, c,Y) = c — Yfic) + (/ - YOJ(X))(X - c) 

where I is the n x n identity matrix, □/(X) is an inclusion function 
for the Jacobian of / evaluated on X, and c is any point in X. 
Note that the vector addition and subtraction and the matrix/vector 
multiplication operations used in K must be computed using interval 
arithmetic. 

We therefore have 

X* C (XP|X(X,C,D) 

for any c G X and nonsingular matrix Y. Thus, KiX, c, Y) can be 
used as an interval Newton operator. Fairly good results can be 
achieved with c = mid(X) and Y = J~\c ) [TOTH85,MITC92], 

In our research, we have found a different method to be superior, 
described in the next section. 


The next two sections discuss two methods for solving these interval 
linear equations. 

4.2.1 Solving the Interval Linear Equation with the Interval 
Inverse 

One method to bound the set of solutions to the interval linear 
equation involves computing the interval inverse. We seek a bound 
on QiM, b): the set of solutions, x, for Mx - b. If M is an n x n 
interval matrix, an interval that bounds QiM, b) is 

Z = M~ l b 


14 In this case, the interval b has a lower bound equal to its upper bounds in each 
coordinate (called a point interval), neglecting inaccuracies in the computation of /. 
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where M~ l is the interval inverse of the interval matrix. Assuming 
M contains no singular matrices, the interval inverse is an interval 
bound on the set 

{m~ l | m G M} 

A simple way of computing the interval inverse is to use the 
interval analog of LU decomposition. That is, we take the LU 
decomposition algorithm [PRES86, pages 31-38] and replace all 
arithmetic operations with their corresponding interval arithmetic 
counterparts (see [MOOR79] for a discussion of interval arithmetic). 
If at any point in the iteration we attempt to divide by an interval 
which contains zero, then we cannot compute the interval inverse, 
and the interval Newton step fails (but see the next section for a way 
to reduce the size of candidate regions without using the interval 
inverse). After enough iterations of the constrained minimization 
algorithm, and assuming the conditions discussed in Section 5 hold, 
candidate regions are usually small enough to make this technique 
effective. 

4.2.2 Solving the Interval Linear Equation with Matrix Itera¬ 
tion 

Another method to bound the set of solutions to the interval linear 
equation involves matrix iteration. The algorithm we present here 
requires an initial bound on the set of solutions; that is it finds a 
bound on the set Q(M, b)f]X where X is a given interval. The 
algorithm is therefore effective at reducing the size of a candidate 
interval in which solutions to an equality constraint are sought, but 
cannot be used to verify solution existence using the theorems in 
Section 5. 15 

Figure 8 contains the code for Linear_Solve, which finds bounds 
on the solution to the linear interval equation. Linear_Solve is 
based on the observation that the i-th equation of the linear system 
Mx = b : 

Mi\X\ + Mi2X2 + * * * + M' in X n — b[ 

implies, for each j such that My ^ 0, that 

bi — MjkXk 

Xi= — m —• 

The interval analog of this equation may therefore be used for each 
interval matrix entry, My, which does not contain 0 to find a bound 
on one of the variables xj. This bound is intersected with the old 
bound on xj yielding an interval which is possibly smaller but no 
larger than it was. Reducing the size of one interval may then further 
reduce the sizes of others as the iteration proceeds. Note that the 
algorithm does not halt when an interval element of M contains 0; 
it just proceeds to the next element which excludes 0. 

An important property of Linear_Solve is that it can be ap¬ 
plied to a nonsquare linear equation, 16 and is therefore useful in the 
“overconstrained” equality constraint for implicit surfaces, and the 
vertex-to-edge and vertex-to-vertex tangency situations of piece- 
wise parametric surfaces (see Appendix A). Linear_Solve can be 
applied in many situations where LU decomposition fails because 
of the singularity of the interval Jacobian matrix. Even when the 
Jacobian matrix is singular at the solution point , Linear_Solve is 
usually effective at reducing the widths of some of the input vari- 


15 This is because, unlike the technique presented in Sections 4.2.1, this technique 
does not bound Q(M, b ) directly, but instead bounds Q(M, b ) f~j X. 

16 That is, the number of equations, m, is unequal to the number of variables, n. 


Linear-Solve(M,fr,x) 

repeat 

loop through rows of M (i = 1,2,..., m) 

loop through columns ofM (J = 1,2, ... ,n) 
if 0 0 My then 

x'j - (bi - Y^. M ikXk)/Mii 

x j < xj n xj 

if xj = 0 return no solution 

endif 

endloop 

endloop 

while there is sufficient improvement inx 


Figure 8: Interval Linear Equation Solution Algorithm: This algorithm 
computes the interval Newton step (first statement of the algorithm in Fig¬ 
ure 6). 

ables. These features are critical in making the singular situations 
described in Section 5 computationally tractable. 17 

The “sufficient improvement” condition mentioned in the algo¬ 
rithm can be implemented as 

w(x +1 ) < a w(x) 

where a typical value of the improvement factor, a, is 0 . 9 . Here x l 
denotes the interval bound computed after i iterations of the repeat 
loop. Specifying a maximum number of repeat iterations also 
limits the amount of computation. 

5 Termination Criteria 

Two theorems in interval analysis specify conditions under which a 
square system of equations contains a unique solution in a region. 18 

Theorem 1 (Krawczyk-Moore Existence) If K(X, c, Y) C X , 

K(X , c , Y) ^ 0, and ||7 — 7m/(X)|| < 1, then there is a unique root 
in X , and pointwise Newton’s method will converge to it. 

Theorem 2 (Linear Interval Equation Existence) If Q(X) c X 
and Q(X ) ^ 0, then there is a unique root in X. 

The conditions implied by these theorems thus lead to acceptance 
criteria, A, for the constrained minimization algorithm. Implemen¬ 
tation of Theorem l’s conditions is clear from the discussion in 
Section 4.1. To verify the conditions of Theorem 2, we use the in¬ 
terval inverse method discussed in Section 4.2.1. We have been able 
to verify solution uniqueness much earlier (i.e., in larger regions) 
using Theorem 2’s test. 

We note the conditions for Theorems 1 and 2 can only be verified 
when the determinant of the Jacobian of the equality constraint 
function, C, is nonzero in some neighborhood of the solution. For 
collision detection, the Jacobian determinant is zero at a solution to 
the contact and tangency constraints in the following situations: 

• the contacting surfaces become tangent but never interpene¬ 
trate. They can even stay tangent for an interval of time. 


17 We prefer the method of matrix iteration described here to a faster method (the 
interval analog of Gauss-Seidel iteration) which involves solving only for the diagonal 
matrix elements after a preconditioning step (see [RATS88]). This method requires 
a square system of equations, and will fail when the Jacobian matrix is singular at 
the solution. Interval computations in the preconditioning step also have the effect of 
increasing the size of the solution set Q(M , b) even before any iteration takes place. 

18 See [TOTH85] for a proof sketch and references for Theorem 1, [SNYD92b] for a 
proof of Theorem 2. 
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T 1« T 2 (t) 


Figure 9: Bounding Sphere Collision Time Bound 

• the surfaces contact instantaneously on a curve or higher¬ 
dimensional region. 

• both — the surfaces contact at an infinite set of points through 
an interval of time. 

In these cases, we can not verify that a unique solution exists and 
must resort to heuristic criteria (the “emergency” criteria of the 
constrained minimization algorithm). That is, we consider surfaces 
to have collided when a bound on Cis sufficiently small in a region . 19 

If a solution is on the boundary of a candidate region, we note that 
the conditions of either theorem will be difficult to verify since the 
set result (e.g., Q(X)) will always slightly spill out of the original 
region X. To solve this problem, we should slightly increase the 
region examined for acceptance, X , via 

X 7 = c + (X — c)(l + 7 ) 

where 7 is a small constant (like .1). We then can perform the 
test on this bigger region X f . Then, even if the solution is on the 
boundary of X , the theorem conditions will eventually be satisfied if 
the Jacobian of C is nonsingular in a neighborhood of the solution. 

6 A Simple Bound for the Time of Collision 

We can save time in the collision algorithm by using a fast algorithm 
to reduce the time interval over which collisions are searched. This 
time bound may also tell us, with a minimum of computation, if 
the two surfaces fail to intersect, obviating the need for further 
computation. The test presented here uses ID minimization (only 
over time) rather than minimization over the 4 or 5 dimensional 
space required to solve the full problem. 

As shown in Figure 9, we compute two bounding spheres around 
each of our parametric surfaces. In the case of rigid motion, this 
may be computed beforehand as a preprocessing step. We specify 
two points for each parametric surface 0\ and O 2 about which to 
compute a bounding sphere. These points should be chosen to 
minimize the size of the resulting bound; using the center of mass 
of the surface is a good choice. The bounding radii are: 

Ri = max ||,vi(mi. vi) - Oi|| 

OlVl) 

R 2 = max \\s 2 (u 2 , V 2 ) — O 21 |- 

( M 2 > v l) 


19 This implies, because of the contact constraint, that the surfaces come within a 
specified constant. 


Note that these bounding radii can be computed using a 2D un¬ 
constrained minimization problem, for which the algorithm of Sec¬ 
tion 3.2 is suitable. 

Using the constrained minimization algorithm, this time on a 
simple ID problem, we then find a bound on the time of collision 
via 

t* = minimum{ t | ||7\(0 + 0% - T 2 (t) - 0 2 \\ < Ri +R 2 } 

*£bo >fil 

t* = - minimumj-f | ||7\(0 + 0 % - T 2 (t) - 0 2 || < R\ + # 2 }- 
t£[to >fil 

We can then replace the [to, h] interval in the full collision mini¬ 
mization problem for that pair of surfaces with [to , t* ], or cull the 
pair of surfaces if no solutions to the ID problem are found. 


7 The Full Collision Algorithm 

The complete algorithm for detecting collisions can now be de¬ 
scribed. The following discussion pertains to a set of parametric 
surfaces; a similar algorithm can be developed for the case of im¬ 
plicit surfaces or parametric/implicit combinations. We are given a 
set of solids defined by a parametric boundary representation, and 
a time interval in which to detect collisions, [to, t\ ]. The following 
steps summarize the final algorithm: 

1. Detect pairs of objects which can possibly collide. For this 
step, we bound each time-varying surface by evaluating an 
inclusion function for its time-varying mapping over [to, t\\. 
More precisely, a bounding box through time on the surface 
S(u, v, t) is given by 

□S([0, 1], [0,1], [to, t{\) 

assuming S is evaluated on the unit square . 20 We can then 
test whether any of the resulting bounding boxes intersect us¬ 
ing highly efficient algorithms from computational geometry 
[SIX82]. All pairs of bounding boxes which do intersect must 
be processed further; the rest are culled. 

2. For rigid bodies, additional object pairs can be culled using the 
bounding sphere test of Section 6 . A variant of this test can 
also be used for deformable surfaces. 

3. If any pairs of objects remain to be processed, we must invoke 
the full constrained minimization algorithm. Here, we distin¬ 
guish between “free” objects and objects already in continuous 
contact, whose contact points are being tracked with the ODE 
solver. For objects already in continuous contact, additional 
constraints are added (Section 3.2) to prevent re-detection of 
the tracked points. All such problems are placed on the initial 
priority queue of the multiple element constrained minimiza¬ 
tion algorithm. 

4. We use local methods, such as Newton’s method, to converge to 
the actual collision point in each solution region which contains 
an isolated collision (i.e., for which the interval existence and 
uniqueness test of Section 5 succeeded). We arbitrarily choose 
the midpoint as the collision point for the rest of the solution 
regions (termed singular solutions in Section 3). 


20 Note that for rigid surfaces we can cache a bounding box on the time-independent 
rigid surface and compute only a bound over time on the resulting rotated and translated 
bounding box. 
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8 Implementing Inclusion Functions 

The collision detection algorithm depends on inclusion functions 
for the time-varying surfaces and their various derivatives. Note 
that the equality constraint for the parametric surface case (Equa¬ 
tion 2) involves derivatives of the time-varying surface mappings 
Si(ui, Vi, t). The interval Newton method then requires an additional 
derivative of the equality constraint with respect to each of the inde¬ 
pendent variables. Interval analysis provides the necessary theory 
for constructing inclusion functions for these functions. 

For simple polynomial surfaces (e.g., bicubic patches or algebraic 
surfaces) interval arithmetic suffices to provide an inclusion func¬ 
tion for the time-independent surface. Toth [TOTH85] has presented 
efficient inclusion functions for Bezier surfaces. Mitchell and Han- 
rahan have proposed a simple stack-based representation of surfaces 
which allows generation of inclusion functions for the surface and 
its derivatives [MITC92]. Inclusion functions for more complicated 
surfaces and their derivatives can also be constructed. We have used 
the system described in [SNYD92a,SNYD92b], which automates 
the construction of inclusion functions (and inclusion functions for 
the derivatives) of any functions formed by the composition of a 
quite powerful set of symbolic operators. 

For physical simulations, the ODE solver computes a represen¬ 
tation of the time behavior of the surfaces. The solver may directly 
compute a continuous representation or it may be later reconstructed 
by point sampling the solver’s results, typically producing a poly¬ 
nomial. Appendix B discusses a method to bound Chebyshev poly¬ 
nomials. 

9 Results 

We have successfully tested this method on a series of collision 
detection examples, including both rigidly moving and deforming 
objects. For example, Figure 12 shows the results of a difficult 
collision detection run in which the contact manifold forms a series 
of disjoint 2D regions. A collection of 59 points was generated 
in the contact region with a simultaneity threshold of 0.001 and 
solution separation distance of 0.04, using 28704 iterations and 
88.81 CPU seconds. 21 While the running time may seem large, the 
problem itself is sufficiently difficult that its running time exceeded 
our threshold of 8 CPU hours without the use of every new technique 
presented in this paper: adding the tangency constraint (rather than 
using the contact constraint alone), sorting by upper bound in the 
constrained minimization algorithm (rather than by lower bound), 
and using Linear_Solve for the interval Newton step (rather than 
the Krawczyk-Moore operator). Figure 1,5, and 12-16 show the 
results of the algorithm for several different time-varying shapes. 

The table in Figure 10 compares running times for a second 
example involving two rotating and translating bumpy parametric 
surfaces which collide at an isolated point. Several solution meth¬ 
ods are compared: LEQN (interval Newton using the linear equation 
solution techniques of Section 4.2), KM (interval Newton using the 
Krawczyk-Moore operator), NIN (without interval Newton), and 
NTAN (without the tangency condition). Since the collision occurs 
at an isolated point, both the LEQ and KM methods were able to ac¬ 
cept a single solution region by verifying the solution existence and 


21 The term iteration refers to an evaluation of the inclusion functions □/ and OF 
(objective function and constraint function) in the constrained minimization algorithm. 
All CPU times are measured on a HP 9000 Series 750 computer. 


| Running Times | 

Example 

Iterations 

CPU (secs) 

LEQN 

6331 

32.67 

KM 

10087 

148.28 

NIN, 7 =le -3 

17395 

8.58 

NIN, 7 =le -4 

29921 

15.46 

NIN, 7 =le -5 

40127 

21.52 

NIN, 7 =le -6 

48187 

23.25 

NIN,NTAN, 7 =le -3 

52307 

14.59 

NIN,NTAN, 7 =le -4 

587711 

169.87 

NIN,NTAN, 7 =le -5 

3822605 

1207.46 


Figure 10: Table of Results for Various Methods: see Section 9 

uniqueness test. The other methods required an accuracy parameter 
for acceptance; we used the simple criterion w(X) < 7 . 

Because we used a prototype system to gather the data, we em¬ 
phasize the importance of iteration count data over CPU time. Our 
system requires the traversal of a complicated data structure for each 
inclusion function evaluation which overwhelms the floating point 
computation actually needed in the function. The interval Newton 
methods are sensitive to this bias, since their implementation re¬ 
quired many symbolic operators. We believe the iteration counts 
shown here to be a reasonable measure of expected running time, if 
the inclusion functions are hand-coded for the surfaces of interest. 

10 Conclusions 

We have presented a robust interval algorithm that can detect col¬ 
lisions between complex curved surfaces. The algorithm handles a 
greater range of situations than previous algorithms. It detects both 
isolated collision points and collision points on contact manifolds. 
It can avoid detection of points close to a set of tracked points with 
specified trajectories. It efficiently handles detection of simultane¬ 
ous collisions between sets of moving objects. The technique is 
practical for simulations involving large numbers of moving and 
deforming objects (see Figures 15 and 16). 

We draw several conclusions from our experimental results. First, 
interval methods, such as [VONH90] and [DUFF92], which do not 
make use of the interval Newton method or the tangency condition 
soon become impractical as we increase the accuracy parameter 
(refer to the NIN,NTAN lines of the table in Figure 10). Interval 
Newton iteration combined with the tangency condition (especially 
using the interval linear equation approach) is very effective at re¬ 
ducing computation. Second, our method can solve the difficult 
problem of detecting collision points on a contact manifold. We 
have found the methods described here to be indispensable, includ¬ 
ing the idea of the tangency constraint, the constrained minimization 
algorithm discussed in Section 3, and the interval linear equation 
approach to interval Newton iteration. 

We note that many areas for improvement remain. Sorting by 
lower bound of the objective function rather than by upper bound 
is more efficient for isolated point collisions. We have noted an 
efficiency gain of a factor of from 1 to 10 in using the lower bound 
for such cases. On the other hand, sorting by lower bound is com¬ 
pletely impractical for detecting collisions on a contact manifold. 
If we know the nature of the collision solution set a priori, we 
can choose the appropriate method. Alternatively, combining the 
two approaches, perhaps by “racing” them in parallel on the same 
problem, may decrease the average running time. We are study¬ 
ing several ways to increase efficiency that involve more optimally 
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choosing the next dimension to subdivide, and determining a sub¬ 
division location other than the midpoint. 
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A Collision Constraints for Piecewise Parametric 
Surfaces 

A piecewise surface is composed of a set of smooth faces, a set of edges where 
these faces meet, and a set of vertices or points where edges meet. Edges 
form the ID boundaries over which the surface is not smooth; vertices are 
the 0D boundaries between smooth edge curves. A data structure containing 
the faces, edges, and vertices of a solid is called its boundary representation . 
For example, the boundary representation of a cylindrical solid contains 
three faces: one cylinder and two circular endcaps, two edges where the 
cylinder and endcap meet, but no vertices. 

To detect collisions between two piecewise surfaces, we must search for 
collisions between each pair of faces, between edges and faces, between 
vertices and faces, etc. The constraints governing collisions are different in 
each of these cases, which we call tangency situations . There are 6 types of 
tangency conditions in a collision between piecewise surfaces as shown in 
Figure 11. Constraints for the face-to-face tangency situation are identical to 
the constraints discussed in Section 2.1. The following paragraphs discuss 
the other tangency situations. 

We must combine all the constrained minimization problems for the 
various possible types of tangency situations. For example, if the surfaces are 
a pair of cylindrical solids, we obtain 25 separate constrained minimization 
subproblems: 9 face-to-face problem, 12 edge-to-face problems, and 4 
edge-to-edge problems. Two tori require only a single face-to-face problem. 
Each problem is then solved simultaneously using the multiple element 
constrained minimization algorithm. 


332 


Edge-to-Face For the edge-to-face case, we have an edge curve, C(s , t), 
which forms a boundary of a surface, S a (u a ,v a ,t), and another surface, 
S(u , v, t). The edge curve is typically formed by evaluating a parametric 
surface S a at a specific value for either the u or v parameter, e.g. 

C(s,t) = S a (s,V UXed ,t) 

where v^ xed is a constant set at one of the extremes of the v interval over 
which S a is evaluated. The contact constraint for edge-to-face collisions is 

C(s , t ) — S(u , v, t) = 0 (8) 

and the tangency constraint by 

^(s,t)-N(u,v,t) = 0 (9) 

os 

where N is the time-varying normal to the surface S. The edge-to-face 
equality constraint can be represented a system of 4 equations in 4 variables. 

To define the incoming collision condition, we need to define what “out¬ 
wardness” means on an edge curve. Assuming all surfaces form the valid 
boundaries of a closed solid, the edge curve C(s , t ) is shared between two 
surfaces S a and Sb- We can therefore define two “outward” directions, given 
by the outward pointing normals to the shared surfaces S a and Sb- For 
example, these outward directions may be defined as 

C°utward- 1 (i 0 = NJs .,,fixed 

c outward- 2 (i0 = Nb (J™ d ,s,t) 

where N a and Nb are the outward normal vectors of the respective surfaces. 

The incoming constraint forces the relative velocity between the surface 
and the edge curve to be in the same direction (using a dot product test) as 
the surface’s normal. The surface’s normal must also face away from at least 
one of the edge curve’s outward directions. The incoming constraint is: 

dS dC 

(—(w,v, t) - (s,t))-N(u,v,t) > 0 and 

m Ot '" 

(~N(u, v, 0 • c outward ' 1 (s, t) > 0 or (10) 

—N(u, v, t) ■ C outward ' 2 (i, () > 0^ . 

Edge-to-Edge The edge-to-edge case involves two edge curves, C\ (s \, t ) 
and C 2 CS 2 , 0- For this case, just a contact constraint is sufficient, given by 
the following system of three equations in three variables: 


system of three equations in a single unknown for the vertex-to-vertex case. 

The incoming collision condition can be derived by defining a number 
of outward directions for the colliding vertex, corresponding to the normal 
vector of each surface containing that vertex. The normal to the surface 
S must face away from at least one of these outward directions, as in the 
edge-to-face case. The relative velocity between the surface and the vertex 
must face in the same direction as the surface’s normal, via 

Ot ot 

Similar systems of inequalities can be derived for situations where a vertex 
collides with an edge or another vertex. 

B Inclusion Functions for Chebyshev Polynomials 


Chebyshev polynomials are a good basis for a continuous representation of 
time behavior. They allow simple control of approximation error, and can be 
differentiated using a simple method to produce a Chebyshev representation 
of the derivative (see [PRES86, pages 158-165] for a discussion of the 
advantages of Chebyshev polynomials, their properties, and algorithms for 
their manipulation). The basis functions for a Chebyshev polynomials are 

T n (x) = cos(n arccos(x)) 

which expand to a series of polynomials of the form 

To(x) = 1 

T\ (x) = x 

T 2 (x) = 2x 2 — 1 

T n +\ (x) = 2xT n (x) - T n _ i(x) n > 1 

The function T n (x) has n + 1 extrema with values of ± 1 at the locations 

7 ri 

Xi = cos(—) i = 0,1 ,..., n 
n 

The i -th extremum of the basis function T n is either a minimum or maximum 
according to the rules 

T ( x / - 1 , if (/ + n) = 1 mod 2 

+i, if (/ + ft) = 0 mod 2 

A Chebyshev approximation of order N is given by specifying N coeffi¬ 
cients Cj, i = 0, 1,... , N — 1, which determine the polynomial 


Ci(si,f)-C 2 (s2,f) = 0. (11) 

To define the incoming collision condition, we define two outward direc¬ 
tions for each edge curve, as in the previous discussion for the edge-to-face 
case. The relative velocity between the edge curves must face in the same 
direction as at least one of the first curve’s outward directions: 

(^(*1,0 - • c° utward " 1 (.s|.f) > 0 or 

%t dt (12) 

(^■(*1,0 - 0 ) • c° utward - 2 (,i, 0 > 0 

Ot ot 

Also, at least one of the outward directions on one curve must face away 
from one of the outward directions of the other curve. A logical combination 
of 6 inequalities is the result. 

Vertex-to-Face, Vertex-to-Edge, Vertex-to-Vertex The vertex-to- 
face case involves a vertex Pit) and a surface S(w,v, t). As in the edge-to-edge 
case, a contact constraint is sufficient, of the form 

P(t) -S(u,v,t) = 0. (13) 


N-l 

C(x) = J2 CiTi ^> + f 
1=1 

Given the order of the Chebyshev approximation function C(x), N, we 
can easily compute an inclusion function for C(x). Let the interval over 
which we are to bound C(x) be given by X = [xo, jci]. As a preprocessing 
step, we first tabulate the locations of the extrema of the basis functions, up 
to some maximum order. (Note that the results can then be used for any 
approximating polynomial.) For each Chebyshev basis function, 7)(x), i = 
0.1..... /V — 1, we first evaluate 7)(xo) and 7)(xi). We then determine 
whether any extrema of 7)(x) occur in [xo, x\] using the tabulated locations 
of the extrema. A lower bound on the basis function over [ xq , xi ], b® is 


f min(7)(xo), 7/(x 1 ), -1), if min of 7)(x) G [x 0 ,xi] 
\ min(T/(xo), Tj(x\ )), otherwise. 


Similarly, an upper bound is 


f max (Ti (xo), 7)(xi), 1), if max of T[ (x) G [xo, xi ] 

\ max(T/(xo), 7/(x 1 )), otherwise. 


where the point P is formed by evaluating a surface at a fixed point in its 
(w, v) parameter space, e.g. 

P(t) = 5 a (« fixed , y fixed 1 1) 

A system of three equations in three unknowns results. Similarly, a system 
of three equations in two unknowns results for the vertex-to-edge case, and a 


The final inclusion function is then 
N-l 

□ C(X)= Y^cAb^b}]*^ 

i= 1 

where operations are computed with interval arithmetic. 
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Scenes from test animations: In the following figure pairs, the upper image is the scene immediately before the collision, while the bottom 
image is the scene at the collision time. Points of contact are shown as white dots, which are uniformly distributed over regions where there are 
line and surface contacts. At the time of collision, surfaces become transparent to make the dots visible. 



Figure 12: Two bumpy objects collide at one Figure 13: A time-varying tube contacts a cush- Figure 14: A wavy object contacts a raised 

point. ion along a curve. checkerboard floor in several flat patches. 


Scenes from “Fruit Tracing”: This animation shows the results of collision detection for a more complicated setting involving hundreds of 
colliding objects. In this animation, moving parametric surfaces representing fruit are collided with a static lobster shape, defined as an implicit 
surface. (Lobster data generated by David Laidlaw, Matthew Avalos, Caltech, and Jose Jimenez, Huntington MRI Center.) 



Figure 15: Colliding dynamic fruits. 


Figure 16: Scene showing lobster shape. 
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